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Abstract 

A steady state mathematical model was developed to 
describe the operation type of multicomponent distillation 
process using a plate-to-plate method involving Murphree 
plate efficiency. The developed model equations were 
applied to a ternary system of acetone-chloroform-methanol 
in a ten-stage distillation column. An algorithm for solving 
the operation type of this ternary distillation problem was 
presented. Numerical solutions were obtained via MATLAB 
using condenser's boundary condition as initial guess. The 
simulated results of the variations of temperature, activity 
coefficient and mole fractions of acetone, chloroform and 
methanol in the liquid, with number of theoretical stages in 
the 10-stage distillation column, conformed to the trends in 
the literature. The convergence behaviour of the method 
employed in this work showed that convergence was 
achieved within 7-10 iterations. 
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Introduction 

Distillation is by far the most widely practiced 
technique to separate mixtures of chemical species in 
petroleum, natural gas, liquid and chemical industries. 
It is a separation process that takes advantage of 
differences in boiling points in a number of 
vaporisation and condensation steps (Solomon, 1981). 
Brooks (1993) defined distillation column as a device 
used to separate components of a mixture using the 
thermodynamic properties of the components. A 
schematic diagram of a typical distillation unit with a 
single feed and two product streams is shown in 
Figure 1. The column itself is adiabatic. The initial 
mixture to be processed, known as the feed, is 
introduced somewhere near the middle of the column 
to a tray known as the feed tray, / . The feed tray 


divides the column into a top (absorption/enriching/ 
rectification) section and a bottom (stripping/ 
exhausting) section. Heat is added to the bottom of the 
column as reboiler heat, Q R , and removed at the top 
as the heat load of the condenser, Q c . Essentially, 
components with lower boiling points will tend 
toward the top of the column as vapour and the 
components with higher boiling points will gravitate 
toward the bottom as liquid. 

Hence, distillate, D , and bottoms, W , are removed 
from the top and bottom of the column respectively. 
Inside the tower, liquids and vapours are always at 
their bubble and dew points respectively so that the 
highest temperatures are at the bottom, the lowest at 
the top. The reboiler can be considered as a theoretical 
stage or tray. Thus, there are internal flows of vapour 
and liquid within the column as well as external flows 
of feeds and product streams, into and out of the 
column. The contact between vapour and liquid on a 
tray engenders heat and mass exchanges conducive to 
the separation goal of a distillation column. The 
entropy and molar flow between the liquid and 
vapour phases are complex phenomena. Convection, 
conduction and heat of vaporisation take place as the 
fluids flow through a stage, exchanging heat and mass. 

Modelling and simulation of distillation columns are 
not new enterprises. The models used for distillation 
columns vary from the simple to the quite complex, 
depending on the nature of assumptions. The models 
described in the literature either contain algebraic 
loops or simplified assumptions that render the models 
ill-equipped for dynamic simulations (Brooks, 1993). 

Multicomponent distillation calculations are classified 
as design or rating methods (Rao et al., 2005). In the 
former, the number of stages has to be determined 
given the recovery of key components whilst in the 
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latter, the composition of products has to be found 
given the number of stages in the enriching and 
stripping sections of the distillation column. Lewis et 
al. (1932) proposed a method of the stage-to-stage 
calculations for multicomponent distillation. To 
determine an approximate distillate composition, they 
considered that components lighter than light key will 
present in significant amounts only in the distillate 
and components heavier than heavy key in the bottom 
product. Assigning arbitrarily small amounts for the 
missing components, they proposed to compute the 
product compositions to initiate the stage-to-stage 
calculation and match the feed stage compositions 
computed from both sections. The discrepancy 
between the feed stages' compositions thus being 
computed was used to revise the products 
compositions. However, this method did not gain 
acceptance and currently there are no simple and 
robust design methods (Fidkowski et al., 1991). To 
avoid the arbitrary assignment of the missing 
components in the products, Thiele et al. (1933) 
proposed a method in which they preassigned the 
number of stages in each section of the column, and 
determined the product compositions. This method 
lends itself well to cast the component and enthalpy 
balances around the stages into elegant tridiagonal or 
block diagonal matrix equation(s). 

Then, the Newton-Raphson method was used to find 
the compositions of products, and liquid and vapour 
streams in the column. The developments in this 
direction have led to the well-known Napthali- 
Sandholm method and its variants. To extend the 
domain of convergence and find multiple solutions, 
physical and mathematical homotopy-continuation 
methods were used. Ketchum (1979) integrated the 
relaxation method, proposed by Rose et al. (1958), into 
a variant of the Napthali-Sandholm method. These can 
be adopted for equilibrium and non-equilibrium stage 
models and reactive distillation. Hass (1992) reviewed 
the developments of these methods. These methods 
are robust and several commercial softwares have 
incorporated them for solving rating problems. The 
above methods are used for design as well, but the 
solutions of several runs of rating problems are 
required. Based on experience, using the Fenske- 
Underwood-Gilliland short-cut method, one obtains 
the near about solution for the design problem. This 
solution can be used to minimize the rating runs 
required to arrive at the solution of design problem. 
This, perhaps, constitutes the widely design method at 
present. 

The design method had been previously revised in the 


context of azeotropic distillation to determine 
rigorously the minimum reflux. Levy et al. (1985) and 
Julka et al. (1993) proposed methods for the direct and 
indirect splits of ternary and multicomponent 
mixtures based on the 'distillation' lines and showed 
that they could be used for any split of multi- 
component systems. Brooks (1993) worked on 
modelling and simulation of a distillation column 
using bond graphs. The structure and the equations 
that represent a sieve-tray distillation column were 
explored using bond graphs. They concluded that 
creating a truly dynamic, rigorous bond graph model 
of a distillation column requires more fundamental 
research. Thong et al. (2000) (2001) proposed methods 
using the 'stage-composition' lines. In a way, these are 
similar to the Lewis-Matheson method, but the match 
is found either geometrically or by a 'tuning' algorithm. 

Jimenez et al. (2002) provided a basis for the 
estimation of liquid-phase activity coefficient of any 
component i , from the molecular volume fractions, O • . 
The liquid-phase activity coefficients can be 
individually differentiated in the combinatorial part, 
which includes the geometric significance for 
combining molecules of different shapes and sizes, 
and the residual part, which includes the energy 
parameters. Rao et al. (2005) worked on simple design 
method for multicomponent distillation columns. 

A simple design method in which the stage-to-stage 
calculations were initiated from the feed stage with the 
feed compositions was presented. The method 
permits the use of component or matrix tray 
efficiencies or the non-equilibrium-stage model to 
estimate the number of stages. The convergence was 
achieved in all the cases within 3-20 iterations taking 
CPU time in the order of 100 th of a second on a PC 
with a Pentium IV processor. Zuzana et al. (2007) 
examined the impact of mathematical model selection 
on prediction of steady state behaviour of a reactive 
distillation column. The objective of their paper was 
to compare the equilibrium (EQ) and non-equilibrium 
(NEQ) models of a reactive distillation column, 
focusing on the phenomenon of multiple steady states. 
It was concluded in their work that the EQ and NEQ 
models showed qualitative different responses to feed 
flow rates disturbances of methanol and butanes, 
predicted different start-up strategy to achieve higher 
steady states. Minh et al. (2009) introduced a 
calculation procedure for modelling and control 
simulation of a condensate distillation column based 
on the energy balance structure. The mathematical 
modelling simulation was accomplished over three 
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phases: the basic nonlinear model, the full-order enriching and stripping sections of the column must 
linearized model and the reduced-order linear model. be equal for all components. 


Recently, Frau et al. (2010) investigated estimation 
model design for a multicomponent distillation column 
with temperature measurements. This approach was 
applied through simulations to an industrial-type six- 
component column. It was concluded that the estimation 
task could be performed with much less (57) ODEs 
compared to the standard three (2849 ODEs) and six 
(17390 ODEs) component Extended Kalman Filters. 
Truong et al. (2011) reviewed some techniques of 
modelling and simulation of distillation columns and 
described a model for a lab-scale binary continuous 
distillation column with the focus was on the dynamic 
behaviour of the product compositions under feed 
disturbances. The simulation results showed that the 
composition responses to disturbances were close to 
the response of a first order system. Popoola et al. 
(2013) did comprehensive review of expert system 
design of crude oil distillation column as an industrial 
application of distillation column using the neural 
networks model. It was observed that the use of 
dynamic mathematical models was challenging due to 
their time dependency; and then they proposed the 
use of back-propagation algorithm involving error 
minimal method to replace convergence problem 
experienced by previous researchers. 

The objective of this study was to develop a steady 
state mathematical model for multicomponent 
distillation column using plate-to-plate method 
involving Murphree plate efficiency and a condition of 
plate matching at the vapour on the feed plate. Hence, 
an algorithm was proposed for solving the operation 
type of multicomponent distillation problem with a 
view to determining the liquid composition and 
temperature profile, variation of activity coefficient 
with number of stages, and the number of iterations 
required to achieve convergence. 

Scope 

Most of the calculation methods used in the studies of 
distillation process were developed by assuming a so- 
called ideal stage. In this study, a steady state 
mathematical model was developed for the operation 
type of multicomponent distillation process by use of a 
plate-to-plate method involving Murphree plate 
efficiency and a condition of plate matching at the 
vapour composition on the feed plate in contrast to the 
0's method proposed by Holland et al. (1983), where 
plate matching implies that the liquid composition on 
the feed plate obtained by calculation steps for both 


The Mathematical Model 


The multicomponent distillation column was modelled 
using the plate-to-plate calculation technique 
involving Murphree plate efficiency and a condition of 
plate matching at the vapour composition on the feed 
plate. Figure 1 shows a typical distillation column, 
whereby the feed F, is introduced on plate / , and 
overhead and bottom products are withdrawn as 
distillate, D, and bottoms, W, respectively. The feed 
tray divides the column into a top (absorption/ 
enriching/rectification) section and a bottom (stripping/ 
exhausting) section. The plates are numbered from the 
top through the feed plate to the bottom of the column, 
as shown in Figure 1. 



FIGURE 1 A TYPICAL DISTILLATION COLUMN. 


The Murphree plate efficiency, Ey^ , is defined by: 

E UV = yj,n-yj,p J = 1;2j M; n = 1,2, , N (1) 


^J,n 


yj,n-yj,v 


where p = n + l for the enriching section and p = N-f 
for the stripping section. 


The operating line of the enriching section is given by: 

y/.n+i = s+l X j ,n + ipl ^ 

where R is the reflux ratio j . Equimolar rates 

of flow of liquid and vapour are respectively assumed 

in the enriching section, that is: 1^ -L 2 - = L n and 

V\=V 2 = = V n+I . 


Using equations (1) and (2) to eliminate y />+ , , and 
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noting that the equilibrium relationship is given by: 

y),n = K j,n X j,n ( 3 ) 

the following equation is obtained: 

yj, n -a 2 x hD (\-E?j) 


X j,n = ' 




so that 


M M 

lL X j,n = Z 

7=1 7=1 


(1 -E%) + E%K m 

y,.n - a 2 x j,p( l - e m ) 

<*,(l ~E") 


. tjMV ^ 

+ h j,n K j,n 


= 1.0 


(4) 


(5) 


where a A = and a 1 = . 

R + l R + l 

Generally, Kj, n is given by a function of both activity 

coefficient (related to composition) and temperature 
thus: 


K M =®(r M ,T„) (6) 


However, K Jn can be regarded as independent of 

composition for an ideal system. The values of 
x j,D ’ and c %2 can be assigned, and yjx is 

determined by the boundary condition at the 
condenser. If the system is ideal, x J n is calculated for 

the given values of E ^ from equation (4) with Kj ^ n , 
which is obtained by solving equation (5). Then, y h2 

is calculated by equation (2), and the stage-to-stage 
calculation is continued toward the feed plate as 
y' = 2,3, , / . If the system is non-ideal, and EfJ is 

given by a complicated function of y* n and y- 

(Yamada et al., 1980), two additional convergence 
loops are considered in the calculations mentioned 
above, where the minor loop is to correct the activity 
coefficient, /j n/ and the major loop is to correct the 

EfJ, as highlighted below. 


The vapour pressures of pure components for a non- 
ideal system are correlated with temperature by 
Antoine equation: 


log p° = A 


B 


(7) 


T + C 

where A , B and C are constants, different for each 
component and T is the temperature in °C. Hence, 

f r» A 


Yj,„ !°g 


-1 


A- 


B 


K j,n ~ 


T + C 


(8) 


The activity coefficient of a ternary system is 
calculated by Margule's equation (Wohl, 1946): 

In Yu = 4,1*2 (! - 2x, ) + 2A x l x x x 2 (1 - ^ ) + A 3 l x 3 2 (l - 2x x ) 


+ 2 A x 3 x x x 3 (l - X] ) - 2A 3 2 x 2 x i ~ 24 3*2*3 + 


(9) 




+ A 2 j + A 2 3 + A 


45,1 


4,2 




(*2*3 “ 2 *1*2*3 ) 


In Yu = 

^3,2 x 3 

( 1 - 

2 x 2 ) + 2A 2 3 

x 2 x 3 (1 - 

x 2 ) ^ A x 2 X\ (l — 

2*2 ) 

+ 2A 2 ,x 

1*2 (! - 

* 2 ) 

2A X ^Xj x^ - 

- 2 A 3 - [ xp 



}(4,2 

+ A, 3 ■ 

■f A 2 

3 + A 3l + A 3 _ 


(xjX 3 - 2 x^X 3 ) 


In r 3 ,n = 

A,3 X \ 

(1- 

2x 3 ) + 2 A 3 j 

* 1*3 (1 _ 

X 3 ) + A 2 3 X 2 ( 1 — 

2 * 3 ) 

+ 2 A 3 2 x 

2 X 3 I 1 - 

-* 3 ) 

| — 2A 2 jX x x 2 

- 24 , 2*1 

X 2 + 



( 10 ) 


( 11 ) 


-(A, 2 + A, 3 


’2,1 


4,3 + 4,i ) _ Q' 


[x x x 2 - 2x x x 2 x 3 ) 


Next, consider the stripping section of the distillation 
column. The operating line of the stripping section is 
given by: 


X j,n+l 


R' Xj, w 

V • “I - 

R r - 1 yj ’ n R' - 1 


( 12 ) 


where R' is the reflux ratio at the bottoms 

The reboiler used in the bottoms is to be an ideal stage, 
E^j =1 . The plate-to-plate calculation in the stripping 

section is continued upward toward the feed plate 
without any complicated procedure, since the values 
of Xj n and y j^- \ to calculate y^ n have already been 

obtained by the previous step. 

The component material balance across the column is 
given by: 

F z j,F = Fx- d + Wxjyy (13) 



Convergence Method for Operation Type 

The so-called operation type is to calculate x j D and 
Xj W for the given operating conditions such as Fxj F 

and R , and number of plates in the enriching and 
stripping sections, n and s respectively. When ideal 
column was assumed, E ^ =1 for all plates and 

components. Holland et al. (1983) proposed the 0's 
method to correct the values of Xj iD and Xj W 

assumed, by satisfying the condition of "plate 
matching" and the component overall material 
balance, where plate matching implies that the liquid 
composition on the feed plate, Xjj , obtained by 

calculation steps for both enriching and stripping 
sections must be equal for all components. However, 
this condition cannot be applied to the plate-to-plate 
method involving Murphree plate efficiency. This is 

because the value of y. + A~y * , determined by 

J ’ \ 7>/ 
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calculation step of the stripping section) cannot be 
obtained in the calculation step of the enriching 
section. Therefore, in the plate-to-plate method 
involving the Murphree plate efficiency, plate 
matching must be made for a condition in which the 
vapour composition leaving the feed plate, y jj , 

calculated from both sections of the column must be 
equal to each other. Thus, equations to obtain the 
corrected values of Xj D and Xj W are as follows: 





X j,D = X j,D(as ) yjJ{W) 

X J’ W \cal L X J,W(as) \[yj,f(D ) 


( 18 ) 


The flow chart of the proposed method is given in 
Figure 2. 


where 6 X and 0 2 are determined by: 



FIGURE 2 FLOW CHART FOR SOLVING THE OPERATION TYPE OF MULTICOMPONENT DISTILLATION PROBLEM. 
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FIGURE 3 TEMPERATURE PROFILE IN THE DISTLLATION COLUMN AT DIFFERENT MURPHREE PLATE EFFICIENCIES 



FIGURE 4 LIQUID-PHASE COMPOSITION OF ACETONE IN ACETONE-CHLOROFORM-METHANOL SYSTEM DOWN THE COLUMN AT 

DIFFERENT MURPHREE PLATE EFFICIENCEIS 



FIGURE 5 LIQUID-PHASE COMPOSITION OF CHLOROFORM IN ACETONE-CHLOROFORM-METHANOL SYSTEM DOWN THE 

COLUMN AT DIFFERENT MURPOHREE PLATE EFFICIENCIES 
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FIGURE 6 LIQUID-PHASE COMPOSITION OF METHANOL IN ACETONE-CHLOROFORM-METHANOL SYSTEM DOWN THE COLUMN 

AT DIFFERENT MURPHREE PLATE EFICIENCIES 



FIGURE 7 VAPOUR-PHASE COMPOSTION OF ACETONE IN ACETONE-CHLOROFORM-METHANOL SYSTEM DOWN THE COLUMN 

AT DIFFERENT MURPHREE PLATE EFFICIENCIES 



FIGURE 8 VAPOUR-PHASE COMPOSITION OF CHLOROFORM IN ACETONE-CHLOROFORM-METHANOL SYSTEM DOWN THE 

COLUMN AT DIFFERENT MURPHREE PLATE EFFICIENCIES 
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FIGURE 9 VAPOUR-PHASE COMPOSITION OF METHANOLL IN ACETONE-CHLOROFORM-METHANOL SYSTEM DOWN THE 

COLUMN AT DIFFERENT MURPHREE PLATE EFFIECIENCIES 


Simulation Result and Discussion 

The model equations developed for the operation type 
of a distillation column were applied to a ternary 
system of acetone-chloroform-methanol, which were 
referred to as components 1, 2 and 3 respectively. The 
model equations were solved using an in-built 
program in MATLAB. An algorithm for solving the 
operation type of multicomponent distillation problem 
is shown in Figure 2. The input data are shown in 
Table 1. 


TABLE 1 INPUT SPECIFICATIONS FOR ACETONE-CHLOROFORM- 
METHANOL SYSTEM. 


Component 

Z J,F 

Antoine's constants 

A. 

B. 

C 

Acetone (1) 

0.2 

7.11714 

1210.595 

229.664 

Chloroform (2) 

0.3 

6.95465 

1170.966 

226.232 

Methanol (3) 

0.5 

8.08097 

1582.271 

239.726 


The constants for calculation of non-ideality are (Perry 
et al., 1999; Prausnitz, 1999): 

4 2 =0.69, 4 3 =0.519, A 2l =0.83, 

4 3 =1.80, 4 a =0.701, 4 2 =0.715 
and Q'=- 0.368. 

In the application of the proposed method, a feed of 
100 mol/h containing 20 mol % methanol, 30 mol % 
chloroform and 50 mol % methanol is to be 
fractionated in a distillation column operating at 
atmospheric pressure of 760 mmHg, with fractional 
recovery of acetone and chloroform in the distillate 


being 0.99 and 0.01 respectively. The distillate is 20 
mol/h. The total number of theoretical stages in the 
column is 10. The parameters used in obtaining 
simulated results are as follows: F= 100 mol/h, D- 20 

mol/h, x Dy =0.99, x D1 =0.01 , N= 10 , 0.4 < EfJ < 1.0 , 
reflux ratios at the top and bottom of the column are 
R=3.0 and R* = 2.5 respectively, and boundary 
condition at the condenser, Y x = [0. 1 5 0.30 0.55] . 

Figure 3 shows the temperature variation in the 
distillation column as one cascade down the column. 
From stages 1 to 2, temperature increases for 
E^J =0.4 and 0.6 and then decreases onwards to 

stage 10 while for E^[ =0.8 and 1.0, a decrease in 

temperature down the column from stages 1 to 10 was 
observed. 

Figures 4 and 7 show the respective increasing liquid- 
phase and vapour-phase compositions of acetone (in 
the ternary system of acetone-chloroform-methanol) 
down the column at different Murphree plate 
efficiencies between 0.4 and 1.0. This is indicative of 
the fact that acetone is the most volatile component. 

The liquid-phase and vapour-phase compositions of 
chloroform and methanol in the ternary system 
decreased down the column at different Murphree 
efficiencies between 0.4 and 1.0 as shown in Figures 5 
and 8; and 6 and 9 respectively. Methanol is the least 
volatile component with chloroform as the 
intermediate. 
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FIGURE 10. VARIATION OF ACTIVITY COEFFICIENT OF ACETONE IN ACETONE-CHLOROFORM-METHANOL SYSTEM WITH STAGES 



FIGURE 11 VARIATION OF ACTIVITY COEFFICIENT OF CHLOROFORM IN ACETONE-CHLOROFORM-METHANOL 

SYSTEM WITH STAGES 



FIGURE 12 VARIATION OF ACTIVITY COEFFICIENT OF METHANOL IN ACETONE-CHLOROFORM-METHANOL 

SYSTEM WITH STAGES 
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Figures 10-12 show the variation of activity coefficient 
of the components with number of theoretical stages 
in the distillation column. It is seen that the activity 
coefficients of acetone and chloroform decrease down 
the column, with the former levelling off for 
0A<E^J <1.0 while the activity coefficient of 

methanol increases down the stage. The activity 
coefficient (non-ideality) explains the escaping 
tendency of a component relative to Raoult's law in 
vapour-liquid equilibrium. Methanol thus shows a 
positive deviation from Raoult's law throughout the 
column. 

In the simulated results presented, convergence was 
achieved within 7-10 iterations, which is indicative of 
the fact that the method of convergence used in this 
study is best appropriate for the ternary system of 
acetone, chloroform and methanol. 

Conclusions 

Multicomponent distillation column has been 
modelled using plate-to-plate method involving the 
Murphree plate efficiency and simulated using 
MATLAB program. The calculations were initiated 
with the boundary condition at the condenser. 
However, further work is required to nullify the 
assumption of constant molar overflow. This removes 
the incorporation of the enthalpy balance equation. 
The method is simple with fast convergence 
achievement but there is a need for integrating a 
search method with the main program for finding an 
optimum feed stage based on some appropriate 
objective function. 

Notation 

A j n Margule's constant (constant of non-ideality), 
dimensionless 

D molar rate of distillate, mol/h 

EfJ Murphree plate efficiency, dimensionless 

F molar rate of feed, mol/h 

K distribution coefficient, dimensionless 

M total number of components, dimensionless 

N total number of theoretical plates in the 

column, dimensionless 

n number of plates in enriching section, 
dimensionless 

W molar rate of bottoms, mol/h 

x t mol fraction of component i in bottoms. 


dimensionless 

y t mol fraction of component i in distillate, 
dimensionless 

z t mol fraction of component i in feed, 

dimensionless 

Greek Symbols 

y t activity coefficient of component i 

6 correction factor for convergence 

Subscripts 

as assumed value 

cal calculated value 

co corrected value 

/ feed plate 

j component 

Superscripts 

stripping section 

s number of plates in stripping section 

saturation condition 
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